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We present the first results from our Post-Newtonian (PN) Smootiied Particle Hydrodynamics 
(SPH) code, which has been used to study the coalescence of binary neutron star (NS) systems. The 
Lagrangian particle-based code incorporates consistently all lowest-order (IPN) relativistic effects, 
as well as gravitational radiation reaction, the lowest-order dissipative term in general relativity. We 
test our code on sequences of single NS models of varying compactness, and we discuss ways to make 
PN simulations more relevant to realistic NS models. We also present a PN SPH relaxation procedure 
for constructing equilibrium models of synchronized binaries, and we use these equilibrium models 
^— ^ ' as initial conditions for our dynamical calculations of binary coalescence. Though unphysical, since 

■ tidal synchronization is not expected in NS binaries, these initial conditions allow us to compare 

' our PN work with previous Newtonian results. 

. We compare calculations with and without IPN effects, for NS with stiff equations of state, 

; I ' modeled as polytropes with F = 3. We find that IPN effects can play a major role in the coalescence, 

accelerating the final inspiral and causing a significant misalignment in the binary just prior to 
, final merging. In addition, the character of the gravitational wave signal is altered dramatically, 

showing strong modulation of the exponentially decaying waveform near the end of the merger. We 
also discuss briefly the implications of our results for models of gamma-ray bursts at cosmological 
distances. 
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I. INTRODUCTION AND MOTIVATION 



On . Coalescing compact binaries with neutron star (NS) or black hole (BH) components provide the most promising 
sources of gravitational waves for detection by the large laser interferometers currently under construction, such as 
O ! LIGO 0, VIRGO II , GEO ||,|, and TAMA |||. In addition to providing a major new confirmation of Einstein's 
O^' theory of general relativity (GR), including the first direct proof of the existence of black holes [^[Ql, the detec- 
^ , tion of gravitational waves from coalescing binaries at cosmological distances could provide accurate independent 

measurements of the Hubble constant and mean density of the Universe |^ . 
^ I Expected rates of NS binary coalescence in the Universe, as well as expected event rates in laser interferometers, 
•rH , have now been calculated by many groups (see for a recent review). Although there is some disparity between 
' various published results, the estimated rates are generally encouraging. Finn & Chernoff calculated that an 
advanced LIGO could observe as many as 20 NS merger events per year. This number corresponds to an assumed 
Galactic NS merger rate Tl ~ 10^^ yr^^ derived from radio pulsar surveys However, later revisions increased 
this number to 72. ~ 3 x 10~^yr~^, using the latest galactic pulsar population model of Ref. This value is 

consistent with the upper limit of 7?, ^ 10~^yr~^ derived (most recently) by Arzoumanian et al. p4| on the basis of 
very general statistical considerations about radio pulsars, and by Kalogera & Lorimer [p^ , who studied constraints 
from supernova explosions in binaries. 

Many calculations of gravitational wave emission from coalescing binaries have focused on the waveforms emitted 
during the last few thousand orbits, as the frequency sweeps upward from ~ 10 Hz to ^ 300 Hz. The waveforms in this 
frequency range, where the sensitivity of ground-based interferometers is highest, can be calculated very accurately by 
performing post-Newtonian (hereafter PN) expansions of the equations of motion for two point masses [p^ . However, 
at the end of the inspiral, when the binary separation becomes comparable to the stellar radii (and the frequency 
is ^ IkHz), hydrodynamics becomes important and the character of the waveforms must change. Special purpose 
narrow-band detectors that can sweep up frequency in real time will be used to try to catch the last ~ 10 cycles of the 
gravitational waves during the final coalescence |17| . These "dual recycling" techniques are being tested right now on 
the German-British interferometer GEO 600 [Q. In this terminal phase of the coalescence, when the two NS merge 
together into a single object, the waveforms contain information not just about the effects of GR, but also about the 
interior structure of a NS and the nuclear equation of state (EOS) at high density. Extracting this information from 
observed waveforms, however, requires detailed theoretical knowledge about all relevant hydrodynamic processes. If 
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the NS merger is followed by the formation of a BH, the corresponding gravitational radiation waveforms will also 
provide direct information on the dynamics of rotating core collapse and the BH "ringdown" (see, e.g., Ref. Q). 

Many theoretical models of gamma-ray bursts (GRBs) have also relied on coalescing compact binaries to provide 
the energy of GRBs at cosmological distances . The close spatial association of some GRB afterglows with faint 
galaxies at high redshifts is not inconsistent with a compact binary origin, in spite of the large recoil velocities acquired 
by compact binaries at birth Currently the most popular models all assume that the coalescence leads to the 

formation of a rapidly rotating Kerr BH surrounded by a torus of debris. Energy can then be extracted either from 
the rotation of the BH or from the material in the torus so that, with sufficient beaming, the gamma-ray fluxes 
observed from even the most distant GRBs can be explained Here also, it is important to understand the 

hydrodynamic processes taking place during the final coalescence before making assumptions about its outcome. For 
example, contrary to widespread belief, it is not clear that the coalescence of two 1.4 Mq NS will form an object that 
must collapse to a BH, and it is not certain either that a significant amount of matter will be ejected during the 
merger and form an outer torus around the central object (see Sec. HI below and Ref. pl[|). 

The final hydrodynamic merger of two compact objects is driven by a combination of relativistic and fluid effects. 
Even in Newtonian gravity, an innermost stable circular orbit (ISCO) is imposed by global hydrodynamic instabilities, 
which can drive a close binary system to rapid coalescence once the tidal interaction between the two stars becomes 
sufficiently strong. The existence of these global instabilities for close binary equilibrium configurations containing 
a compressible fiuid, and their particular importance for binary NS systems, was demonstrated for the first time by 
Rasio & Shapiro (Ref. |^^, hereafter RSl-3 or collectively RS) using numerical hydrodynamic calculations. These 
instabilities can also be studied using analytic methods. The classical analytic work for close binaries containing 
an incompressible fluid (e.g., Ref. was extended to compressible fluids in the work of Lai, Rasio, & Shapiro 

(Ref. [Q, hereafter LRSl-5 or collectively LRS). This analytic study conflrmed the existence of dynamical instabilities 
for sufficiently close binaries. Although these simplified analytic studies can give much physical insight into difficult 
questions of global fluid instabilities, 3D numerical calculations remain essential for establishing the stability limits 
of close binaries accurately and for following the nonlinear evolution of unstable systems all the way to complete 
coalescence. 

A number of different groups have now performed such calculations, using a variety of numerical methods and 
focusing on different aspects of the problem. Nakamura and collaborators (see ||2^] and references therein) were 
the first to perform 3D hydrodynamic calculations of binary NS coalescence, using a traditional Eulerian finite- 
difference code. Instead, RS used the Lagrangian method SPH (Smoothed Particle Hydrodynamics). They focused 
on determining the ISCO for initial binary models in strict hydrostatic equilibrium and calculating the emission of 
gravitational waves from the coalescence of unstable binaries. Many of the results of RS were later independently 
confirmed by New & Tohline and Swesty et al. ||2^, who used completely different numerical methods but also 
focused on stability questions, and by Zhuge, Centrella, & McMillan ||2^,|2^, who also used SPH. Davies et al. and 
Ruffert et al. |^,^ have incorporated a treatment of the nuclear physics in their hydrodynamic calculations (done 
using SPH and PPM codes, respectively), motivated by cosmological models of GRBs. All these calculations were 
performed in Newtonian gravity, with some of the more recent studies adding an approximate treatment of energy 
and angular momentum dissipation through the gravitational radiation reaction [p3[ . 



Zhuge et al. |2^,|29[ and Ruffert et al. |3l|j3^ also explored in detail the dependence of the gravitational wave signals 
on the initial NS spins. Because the viscous timescales for material in the NS is much longer than the dynamical 
timescale during inspiral, it is generally assumed that NS binaries will be non-synchronized during mergers. It is 
generally found than non-synchronized binaries yield less mass loss from the system, but very similar gravity wave 
signals, especially during the merger itself when the gravity wave luminosity is highest ||3ll| . 

All recent hydrodynamic calculations agree on the basic qualitative picture that emerges for the final coalescence. 
As the ISCO is approached, the secular orbital decay driven by gravitational wave emission is dramatically accelerated 
(see also LRS2, LRS3). The two stars then plunge rapidly toward each other, and merge together into a single object 
in just a few rotation periods. In the corotating frame of the binary, the relative radial velocity of the two stars always 
remains very subsonic, so that the evolution is nearly adiabatic. This is in sharp contrast to the case of a head-on 
collision between two stars on a free-fall, radial orbit, where shock heating is very important for the dynamics (see, 
e.g., RSI and Ref. @). Here the stars are constantly being held back by a (slowly receding) centrifugal barrier, and 
the merging, although dynamical, is much more gentle. After typically 1 — 2 orbital periods following first contact, 
the innermost cores of the two stars have merged and the system resembles a single, very elongated ellipsoid. At this 
point a secondary instability occurs: mass shedding sets in rather abruptly. Material (typically 10% of the total 
mass) is ejected through the outer Lagrange points of the effective potential and spirals out rapidly. In the final stage, 
the inner spiral arms widen and merge together, forming a nearly axisymmetric torus around the inner, maximally 
rotating dense core. 

In GR, strong-field gravity between the masses in a binary system is alone sufficient to drive a close circular orbit 
unstable. In close NS binaries, GR effects combine nonlinearly with Newtonian tidal effects so that the ISCO should 
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be encountered at larger binary separation and lower orbital frequency than predicted by Newtonian hydrodynamics 
alone, or GR alone for two point masses. The combined effects of relativity and hydrodynamics on the stability of 
close compact binaries have only very recently begun to be studied, using both analytic approximations (basically, PN 
generalizations of LRS; see, e. g., ||35| , |3^ , as well as numerical calculations in 3D incorporating simplified treatments 
of relativistic effects (e.g., 38 ). Several groups have been working on a fully relativistic calculation of the final 
coalescence, combining the techniques of numerical relativity and numerical hydrodynamics in 3D [ ^9|j40[| . However 
this work is still in its infancy, and only preliminary results of test calculations have been reported so far. It should 
be noted that IPN calculations performed by Taniguchi and collaborators to study the location of the ISCO for 
corotating and irrotational binaries find that the ISCO moves inwards as post-Newtonian corrections are increased, 
due primarily to the effect of IPN potentials with momentum-based source terms present in the system. Similarly, 
Buonanno and Damour find that the ISCO for point masses in a binary under GR moves inwards with increasingly 
massive objects. 

The middle ground between Newtonian and fully relativistic calculations is the study of the hydrodynamics in PN 
gravity. Formalisms exist describing not only all lowest-order corrections (IPN) to Newtonian gravity, but also the 
lowest-order (2.5PN) effects of the gravitational radiation reaction Such calculations have been undertaken by 

Shibata, Oohara & Nakamura using an Eulerian grid-based code, and more recently by Ayal et al. using SPH. 
This paper is the first of a series in which we will present a comprehensive study of the hydrodynamics of compact 
binary coalescence using a new PN version of a parallel SPH code which we have been developing over the past two 
years. This work will be the natural extension to PN gravity of the original Newtonian study by RS. 

PN calculations of NS binary coalescence are particularly relevant for stiff NS EOS. Indeed, for most recent stiff 
EOS, the compactness parameter for a typical 1.4 Mq NS is in the range GM/Rc^ ~ 0.1 — 0.2, justifying a PN 
treatment. After complete merger, an object close to the maximum stable mass is formed, with GM/Rc^ ~ 0.3 — 0.5, 
and relativistic effects become much more important. However, even then, a PN treatment can remain qualitatively 
accura te if th e final merged configuration is stable to gravitational collapse on a dynamical timescale (see the discussion 



in Sec. HID). Most recent theoretical calculations (e.g., the latest version of the Argonne/Urbana EOS; see Ref. [[17[ 



and a number of recent observations (e.g., of cooling NS; see Ref. [^) provide strong support for a stiff NS EOS. In 
this paper we represent NS with stiff EOS by simple polytropes with an adiabatic exponent P = 3 (i.e., the EOS is of 
the form = fcrf , where is the pressure, is the rest-mass density, and fc is a constant; see RS2 and LRS3, who 
obtain F ~ 3 for the best polytropic fit to recent stiff NS EOS). 

The most significant problem facing PN hydrodynamic simulations is the requirement that all IPN quantities be 
small compared to unity. Unfortunately, this precludes the use of realistic NS models. Shibata, Oohara & Nakamura 
1^ computed IPN mergers of polytropes with F = 5/3 and a compactness GM/Rc^ = 0.03, leaving out the effects of 
the gravitational radiation reaction. Ayal et al. performed calculations for polytropes with F = 1.6 or F = 2.6 and 
compactness values in the range GM/Rc^ ~ 0.02 — 0.04, including the effects of the gravitational radiation reaction. 
For comparison, a realistic NS of mass M = IAMq and radius R = 10 km has GM/Rc^ — 0.2, i.e., about an order of 
magnitude larger. Unfortunately, performing calculations with artificially small values of GM/Rc^ also has the side 
effect of dramatically inhibiting the radiation reaction, which scales as {GM/ Rc^)^-^ . 

Our PN SPH code combines a new parallel version of the Newtonian SPH code used by RS with a treatment of 
PN gravity based on the formalism of Blanchet, Damour, and Schafer (Ref. hereafter BDS). Our calculations 
include all IPN effects, as well as a PN treatment of the gravitational radiation reaction. We have also developed 
a relaxation technique by which accurate quasi-equilibrium configurations can be calculated for close binaries in PN 
gravity. These serve as initial conditions for our hydrodynamic coalescence calculations. In addition, we present in 
this paper a simple solution to the problem of suppressed radiation reaction for models of NS with unrealistically low 
values of GM/Rc^. 

The outline of our paper is as follows. Section II presents our numerical methods, including the description of 
our new PN SPH code, a discussion of the advantages of using SPH for this work, and the steps taken to make our 
results as realistic as possible. More details on the methods and our treatment of the initial conditions are given in 
the appendices. Section HI presents our initial results, based on two large-scale simulations, with and without IPN 
effects. These are performed for synchronized initial binaries containing two identical polytropes with F = 3. In future 
papers, we will study systematically the dependence of these results on the NS EOS (by varying GM/Rc^ and F, and 
using more realistic, tabulated EOS for nuclear matter at high density), the NS spins (allowing for nonsynchronized 
initial conditions), and the binary mass ratio. Motivation for this future work as well as a brief summary of our 
present results are presented in Section IV. 



II. NUMERICAL METHOD 
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A. The SPH Code 



Smoothed Particle Hydrodynamics (SPH) is a Lagrangian method ideally suited to calculations involving self- 
gravitating fluids moving freely in 3D. The key idea of SPH is to calculate pressure gradient forces by kernel estimation, 
directly from the particle positions, rather than by finite differencing on a grid (see, e.g., ]49|] , for recent reviews on 
the method). SPH was introduced more than 20 years ago by Lucy, Monaghan, and collaborators who used 
it to study dynamical fission instabilities in rapidly rotating stars. Since then, a wide variety of astrophysical fluid 
dynamics problems have been tackled using SPH (see Ref. |5l| ] for an overview). 

Because of its Lagrangian nature, SPH presents some clear advantages over more traditional grid-based methods 
for calculations of stellar interactions. Most importantly, fluid advection, even for stars with a sharply defined surface 
such as NS, is accomplished without difficulty in SPH, since the particles simply follow their trajectories in the fiow. In 
contrast, to track accurately the orbital motion of two stars across a large 3D grid can be quite tricky, and the stellar 
surfaces then require a special treatment (to avoid "bleeding"). SPH is also very computationally efficient, since it 
concentrates the numerical elements (particles) where the fluid is at all times, not wasting any resources on emty regions 
of space. For this reason, with given computational resources, SPH provides higher averaged spatial resolution than 
grid-based calculations, although Godunov-type schemes such as PPM typically provide better resolution of shock 
fronts (this is certainly not a decisive advantage for binary coalescence calculations, where no strong shocks ever 
develop). SPH also makes it easy to track the hydrodynamic ejection of matter to large distances from the central 
dense regions. Sophisticated nested-grid algorithms are necessary to accomplish the same with grid-based methods. 

Our simulations were performed using a modified version of an SPH code that was originally designed to perform 
3D Newtonian calculations of stellar interactions (see Ref. and RSI). Although the fluid description is completely 
Lagrangian, the gravitational field in our code (including PN terms) is calculated on a 3D grid using an FFT-based 
Poisson solver. Our Poisson solver is based on the FFTW of Frigo & Johnson |^, which features fully parallelized 
real-to-complex transforms. Boundary conditions are handled by zero-padding all grids, which has been found to 
produce accurate results and to be the most computationally efficient method (see also [p7|). Since we are primarily 
interested in the gravitational wave emission, which originates mainly from the inner dense regions of NS mergers, 
we fix our grid boundaries to be ±4 NS radii in all directions from the center of mass. Particles that fall outside 
these boundaries are treated by including a simple monopole gravitational interaction with the matter on the grid. 
Our code has been developed on the SGI/Cray Origin 2000 parallel supercomputer at NCSA. MPI (the Message 
Passing Interface) reduces the memory overhead of the code by splitting all large grids among the processors. All 
hydrodynamic loops over SPH particles and their neighbors have also been fully parallelized using MPI, making our 
entire code easily portable to other parallel supercomputers. The parallelization provides nearly linear speedup with 
increasing number of processors up to ^ 10, with a progressive degradation for larger numbers. 

For more details on the Newtonian version of the code, and extensive results from test calculations, see Ref. fs^ . 



To investigate the hydrodynamics of NS binary coalescence beyond the Newtonian regime, the equations of RS 
were modified to account for PN effects described by the formalism of BDS, converted into a Lagrangian, rather than 
Eulerian form. The main equations and definitions of quantities appearing in the BDS formalism are summarized 
briefly in Appendix A. The formalism is correct to first (IPN) order, with all new forces calculated from eight additional 
Poissson-type equations with compact support, allowing for the computation of all IPN terms using the same FFT- 
based convolution algorithm as for the Newtonian Poisson solver. PN corrections to hydrodynamic quantities are 
calculated by the SPH method, i.e., by summations over particles. Dissipation of energy and angular momentum by 
gravitational radiation reaction is included to lowest (2.5PN) order, requiring the solution of one additional Poisson- 
type equation. 

The key changes to the BDS formalism involve a conversion to quantities based on SPH particle positions, rather 
than grid points. In the discussion that follows, a and b refer to quantities defined for particles, or particle neighbors, 
while i and j are spatial indices. The rest-mass density is calculated at each particle position as a weighted sum over 
the masses of neighboring particles, 



where TOq is the rest mass of particle a, and the weights are given in terms of a smoothing kernel W{f, h) by 



B. The Blanchet, Damour, and Schafer PN Formalism with SPH 




(1) 
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Here ha is the smoothing length for particle a, which is updated after every iteration so as to keep the number 
of neighbors as close as possible to a designated optimal value N^. The form of the SPH kernel W used in our 
calculations is the same standard third-order spline used by RS (and most other current implementations of SPH). 
This kernel is spherical and goes to zero for r > 2h. 

The total mass-energy of the system can be calculated as 

Mg^ f d''xr4l + S) = Y,mb{l + Sb), (3) 



where 6 is a IPN correction defined in Eq. (A12), while the total rest mass is 

M = / (fx n = ^ rub. (4) 



b 

In our simple polytropic models of NS, the pressure is calculated from the local density as 

pi-^ ^ kairi^yf (5) 

where ka is a function of the specific entropy of the particle, and F is the adiabatic exponent. The standard Newtonian 
pressure force is given by 

-E'-'^^ + ^jv.W.. (6) 

In the absence of artificial viscosity ( AV) , entropy is conserved and ka is constant throughout the simulation. If AV 
is included, the first law of thermodynamics can be expressed as 



dka r - 1 

^ b 



with a corresponding force on particles given by 

F^^"" = -Y,mbn(^a,b)'^,Wab- (8) 

b 

The expression for ^i^a,b) depends on the particular choice of AV. We use the form proposed by Balsara which 
gives 

/ pi-) p(b) \ 

^ [-^2 + ^ ) + P'^'l.b)). (9) 

where ^(a,b) is a measure of the rate of convergence in the flow. The exact definition of /i(a.6) can be found in Eqs. (15- 
19) of Lombardi et al. [Q, who also show that the optimal choice of the numerical coefficients is a' = [3' = T /2. 
This form was shown to handle shocks properly and minimize the amount of spurious mixing and numerical shear 

viscosity. 

For the FN pressure force (Eq. A17), we now find 

ppress ^ ^ 2_)i^F^Vdro ^ pAV ^ _ (10) 

C C 

where a/c^ <C 1 is a PN correction. In the calculation of dia = (2 — 3T)diU^ — ^diw"^, we must take a derivative of 
the local dynamic velocity-squared field, which we do by SPH summations, i.e., we first write 

d,iw^) = -{d^ir^w^) - w^a.n), (11) 

and we then calculate the derivative terms as 
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dA"^ ^Y.'^bd.Wab, (12) 
b 

ft(r,u.2)(") = J2 mbwff,)d^Wab. (13) 

b 

The nine Poisson-type equations in the fuU PN formahsm of BDS are aU solved by the same FFT convolution 
method. All 3D grids used by the Poisson solver are distributed among the processors in the z-direction. Real-to- 
complex transforms are computed using the rfftwnd_mpi package of the FFTW library [Esl. The source terms 



of the Poisson equations that do not contain density derivatives, Eqs. ( A5 , A6|, Aq) , are laid down on the grid by a 



cloud-in-cell method. All integrals over the density distribution are converted into sums over particles, e.g., 

Source terms containing density derivatives are calculated by finite differencing on the grid, rather than by SPH-based 
derivatives at particle positions. This has two benefits. First, for integrals of the type 

j3 



$ = J d-'x d,r, . . . (15) 

we cannot convert directly from a volume integral to a sum over discrete particle masses. Second, it guarantees 
that the volume integral of the source term vanishes in Eqs. ( A7, A16| ), as it should. Derivatives of the potentials 



are computed by finite differencing on the grid, and then interpolated between grid points to assign values at SPH 
particle positions. 

The calculation of the quadrupole tensor and its derivatives are unchanged in a Lagrangian formulation. When 
calculating the third derivative, however, we found it advantageous to take the SPH expression for the second derivative 
of the quadrupole tensor (RSI), 

Q^J = E M^l'^^f + xf^djU, + xfd.U,), (16) 

b 

and numerically differentiate once with respect to time. The resulting expression differs from the third derivative 
expression given in BDS by a term of O (f^/c^), but all radiation reaction terms into which it enters already contain 
factors of O (w^/c^). While only approximate, this method proved more stable since it does not require the numerical 
evaluation of several second derivatives on a grid. 

For calculations in which we include the radiation reaction, but ignore IPN corrections, all terms containing a factor 
of 1/c^ in Appendix A can be ignored. In th is cas e our equations reduce to those of the purely Newtonian case, with 



two exceptions. First, we include F['^°-'^ (Eq. A19) in the SPH equations of motion, replacing Eq. (A21) by 

w, = -^ + m, + Fr'''. (17) 

Second, the relationship between the particle velocity v and momentum w is given by 

4 G [■^1 

v, = w, + --QfJw,. (18) 
This has beeen shown [js^ to give the correct energy loss rate as predicted by the classical quadrupole formula, 

f4§<«S'og'>- as) 

Ignoring IPN terms reduces the number of Poisson equations to be solved per iteration from nine to two. The obvious 
advantage is a proper handling of the dissipative PN effects, while leaving the hydrodynamic equations in a simple 
form that can be directly compared to the Newtonian case. In addition, because the corrections are 0{v^ /c'), the 
radiation reaction terms always remain small, even when IPN corrections would be large. 

We have performed a number of test calculations to establish the accuracy of our treatment of PN effects in the SPH 
code. These include tests for single rotating and nonrotating polytropes in PN gravity, which we have compared to 
well-known analytic and semi-analytic results ||56[|. In particular, we have verified that our code repr oduces correctly 



the dynamical stability limit to radial collapse for a single PN polytrope with P = 5/3 (see Sec. B 1). 
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C. A Hybrid 1PN/2.5PN Post-Newtonian Formalism 



Throughout this paper, unless otherwise specified, we use units in which Newton's gravitational constant G, and 
the rest mass M and radius i? of a single, spherical NS are set equal to unity. In Newtonian physics, this leads to a 
scale- free calculation (RS). When we include PN effects, specifying the physical mass and radius of the NS then sets 
the value of the speed of light c, and the magnitude of all PN terms. In our units, the compactness ratio GM/Rc^ of 
a NS is expressed simply as 

The equations of BDS assume that all IPN corrections are small. As mentioned in Sec. 1, this places a rather severe 
constraint on the allowed NS mass and radius, since 

1 / M \ /I5km\ , , 



If, for example, we estimate the potential at the center of the star as J7,/c^ ~ 1.5/c^ — 0.2 1 (Eq. A5), which is 
appropriate for F = 3 models, we find that our "first-order" correction term a/c^ (Eq. AlO| ), with F = 3 and no 
internal motions, is 

4-(2-3F)%--7%~-1.5. (21) 

This is clearly problematic since the derivation of the BDS formalism assumes that |a|/c^ ^ 1. For a fixed radius of 
15 km and F = 3, a NS mass < 0.9 A/q, or 1/c^ < 0.09 is required to keep |a| < 1. This problem is less severe for a 
lower value of F, since the coefficient of a is then smaller. For F — 5/3, we have a — — 3C/*, but these configurations are 
kno wn t o be unstable against gravitational collapse for compactness parameters 1/c^ ^ 0.14 (See, e.g., Ref. |5^, and 



Sec. Bl). These problems are the reason why previous PN hydrodynamic simulations of NS binary coalescence have 
used unrealistic NS models with low masses and large radii. In practice, we find that we cannot calculate reliably NS 
mergers including IPN corrections, unless 1/c^ < 0.05, or c > 4. With such a small compactness parameter, radiation 
reaction effects would then be suppressed by a factor ~ 2^ = 32. 

Recognizing that the IPN and 2.5PN terms describe essentially independent phenomena, and that the proper 
form for energy and angular momentum loss holds even if IPN corrections are ignored, we adopt a hybrid scheme. 
Specifically, in this paper, we set c = 4.47 = cipN for all IPN corrections, which is unphysically large, but we use 
a physically realistic value of c = 2.5 = C2.5PN for the 2.5PN corrections, corresponding, for example, to a NS mass 
M = 1.5 Mq with radius R = 13.9 km. We feel that this hybrid formulation provides a reasonable trade-off between 
physical reality and the limitations of the IPN approximation. 

Note that this method should better extrapolate toward physical reality, compared with unrealistically undercom- 
pact NS models. If -I- 2.5PN simulations are interpreted as taking the limit cipn 00 for the IPN corrections, we 
see that by reducing the compactness in both the -I- 2.5PN and -I- 1 -I- 2.5PN cases, the value of C2.5PN is fixed at 
an unphysical value while cipN is varied, which can never truly extrapolate to the physical case. By setting C2.5PN 
to a realistic physical value while varying cip^v, we may be able to extrapolate our results toward a correct physical 
limit. However, a disadvantage of this approach is that it does not allow for direct quantitative comparison with 
full GR simulations of binary NS coalescence. In these simulations, which essentially handle corrections to all orders 
simultaneously, separation into various PN orders has no meaning. 



D. Initial Conditions 



In addition to its normal use for dynamical calculations, our SPH code can also be used to construct hydrostatic 
equilibrium configurations in 3D, which provide accurate initial conditions for binary coalescence calculations. This 
is done by adding artificial friction terms to the fluid equations of motion and forcing the system to relax to a 
minimum-energy state under appropriate constraints (RS). The great advantage of using SPH itself for setting up 
equilibrium solutions is that the dynamical stability of these solutions can then be tested immediately by using them 
as initial conditions for dynamical SPH calculations. Very accurate 3D equilibrium solutions can be constructed using 
such relaxation techniques, with the virial theorem satisfied to better than 1 part in 10'^ and excellent agreement 
found with known quasi-analytic solutions in both Newtonian (LRSl, LRS4, RS2) and PN gravity |Q. The careful 
construction of accurate quasi-equilibrium initial conditions is a distinguishing feature of both our previous Newtonian 
calculations (RS) and our new PN calculations of binary coalescence. In contrast, most other studies have used very 
crude initial conditions, placing two spherical stars in a close binary orbit, and, for calculations that went beyond 
Newtonian gravity, adding the inward radial velocity for the inspiral of two point masses. As demonstrated in 
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Sec. [II A, this leads to a significantly slower inspiral rate. Moreover, spurious fluid motions are created as the stars 
respond dynamically to the sudden appearance of the strong tidal force. These can in turn corrupt the gravitational 
radiation waveforms. Spurious velocities have additional effects in the full IPN case, where spuriou s motions ente r 
repeatedly into the evolution equations, by propagating through the IPN quantities a, f3, and S (Eqs. AlC , All , A12 ) . 
A specific cause of worry is the influence of velocities adding to S, which affects not only the self-gravity of the stars, 
but also their mutual gravitational attraction. 

We have developed a method, described in detail in Appendix B, that allows for more realistic initial conditions 
for FN synchronized binaries. It reduces dramatically the initial oscillations around equilibrium when the dynamical 
calculation is started. Since the method varies rather significantly between the full IPA'' + 2.5PN formalism and the 
Newtonian with radiation reaction formalism, we handle the two cases separately. For the Newtonian case, radiation 
reactio n pla ys no role in the relaxation, entering only into the initial values of the velocity and momentum, v and 
w (Eq. A22| ), upon the start of the dynamical run. The IFN case is considerably more complicated, requiring the 
construction of static single star models, which are then input into a FN binary relaxation scheme. 



III. RESULTS 



We have performed two large-scale hydrodynamic simulations of NS binary coalescence, with and without the 
IFN correction terms of BDS. Both simulations included radiation reaction throughout the entire run, treated in the 
formalism of Appendix A. Hereafter, we refer to these two runs as the Newtonian (N) and post-Newtonian (FN) runs, 
noting that the N simulation did include 2.5FN effects. 

For both runs, we used 50,000 particles per NS (total of 10^), with a F = 3 polytropic EOS. The two NS are 
identical. Synchronized rotation was assumed in the initial condition. The optimal number of neighbors for each SFH 
particle was set to 100. Shock heating, which plays a completely negligible role in the case studied here, was ignored 
and therefore the SFH AV was turned off. All Foisson equations were solved on grids of size 256'^, including the added 
space necessary for zero-padding. For the IFN run, we used a compactness parameter = 0.05 (see Sec. IIC). 

In both runs, we used C2.5PN = 2.5 in calculating radiation reaction terms. The N run required a total of 600 CFU 
hours and the FN run 1200 hours on the NCSA Origin2000, including the relaxation phase. Farticle plots illustrating 
qualitatively the evolution of the system are shown in Fig. |l| (N run) and Fig. |^ (FN run). 



A. Dynamical Instability and the Inspiral Process 

It was shown by RS and LRS that equilibrium configurations for close binary NS become dynamically unstable 
when the separation r is less than a critical value. For Newtonian, synchronized, equal-mass binaries with F = 3, the 
ISCO is at r = 2.95 R. Furely Newtonian calculations for binaries starting from equilibrium configurations with a 
separation larger than this value will show no evolution in the system. Binaries starting from a smaller separation, 
though, are dynamically unstable, and coalesce within a few orbital periods, even without the energy and angular 
momentum loss due to radiation reaction (RS, [p6| , p7| ). 

In simulations with radiation reaction included, coalescence will always be the end result. The limiting factor on 
how large to make the initial separation is the computing time required for the binary orbit to slowly spiral inward. 
Ideally, one should make sure that the stars are in quasi-equilibrium when the orbit approaches the ISCO and the 
inspiral timescale undergoes a shift from the slow radiation-reaction timescale to the much faster dynamical timescale. 

Since the effective gravitational attraction between two stars is increased by FN effects, we expect the ISCO to move 
outwards when IFN corrections are included. This was demonstrated by Lombardi, Rasio, & Shapiro who used 
the same energy variational method as LRS to find equilibrium configurations for binary NS models including IFN 
corrections. Taking into account these results, we used an initial separation of rp = 3.1 R for our N simulation, and 
r — 4.0 R for the FN simulation. As a consequence, there is an ambiguity in the relative time between the two runs, 
which we resolve by adjusting the initial time of the N run such that the maximum gravity wave luminosity occurs 
at the same time in both the N and FN runs. This was found to require shifting the time in the N run backwards so 
that it starts sX t = —13, while the FN run starts a.tt — Q. 

In Fig. H, we show the evolution of the center-of-mass binary separation during the initial inspiral phase for our N 
and FN runs. Fig. ^ shows the inspiral phase of the N run, as well as the inspiral tracks predicted by the classical 
quadrupole formula for two point masses, and by the methods of LRS3 for two corotating spheres and two ellipsoids. 
We note that the results of LRS3 predict for extended objects a significantly more rapid inspiral rate, which is 
confirmed by the numerical run. In addition, we note that the approach of the ISCO is clearly visible in the plots. 
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where the inspiral rate switches from the slow radiation-reaction-driven orbital decay to the faster dynamical infall. 
This appears to happen at r ~ 2.7 R in the N case, in good agreement with previous results. 

Comparing the PN run to the N run, we see that the stability limit must lie at a considerably larger separation. 
This agrees with the results of Lombardi et al. [Q, who find that PN corrections not only move the ISCO outward, 
but also flatten out the equilibrium binary energy curve E{r) near the stability limit (where E{r) reaches a minimum). 
Following the arguments of LRS3, we conclude that unstable inspiral begins when the differential change in binary 
energy as a function of separation becomes smaller than the energy loss rate to gravitational radiation. The condition 
for unstable inspiral, expressed as 

dE_ dEgw ( dr\ 

dr ^ dt \dt) ' ^ ' 

is then encountered further outside the ISCO (as determined for binaries in strict equilibrium), since PN corrections 
decrease the left-hand side. This effect can also be seen in the results of Ayal et al. by careful examination of 
their Fig. 5a. Even though the binary separation in their PN run has a large initial oscillation, caused by the use of 
non-equilibrium initial conditions, it still converges at a much more rapid rate than in their corresponding Newtonian 
model. 

Even though the effective stability limits of N and PN binaries differ by a large amount, their actual inspiral 
velocities are essentially the same from the moment of first contact, at a separation of r ~ 2.5 i?, until the merger 
of the NS cores. The only significant difference is the break in the inspiral velocity for the N run at t ~ 20, which 
occurs as the cores start to come into direct contact with each other. The lack of this feature in the PN run will be 



explained in Sec. [II B 



B. Coalescence 

In Fig. ^, we show the time evolution of the maximum density in both runs. The maximum density is at the 
center of either star initially, but it shifts eventually to the center of the merger remnant. The initial oscillations with 
a period of T ~ 2 — 3 correspond to the fundamental radial pulsations of the polytropes, and represent the errors 
resulting from small departures from strict equilibrium in the initial conditions. We see that Sp/ p ~ 0.01 and 0.05, 
respectively, for the N and PN runs, which provides a measure of the numerical accuracy of the initial conditions. 

As the binary system contracts to separations of r ^ 2.7 R, we see a rather sudden and rapid decrease in the 
maximum density found at the core of each star, corresponding closely with the moment of first contact of the two 
stars, after which the cores get tidally stretched. For the PN run, this follows a gradual increase in the average 
density maximum, which is caused by the contraction of each NS in response to the growing gravitational potential 
of its companion, rather than a p ure tidal effect. This effect, which seems to result primarily from the weakening 



of th e pre ssure force in Eq. ( A17 ) as a becomes more negative in response to the growing gravitational potential 
(Eq. AlO), was also seen by Ayal et al. in one of their runs (Ref. [Q, see their Fig. 6, run P3). When the center of 
mass separation reaches a value of r ~ 2.0 i? the maximum density stops decreasing, turning around and increasing 
sharply as the cores come into direct contact and merge. 

In Fig. ^ we show the gravity wave signatures of both runs. The waveforms in the two polarizations of gravitational 
radiation are calculated for an observer at a distance d along the rotation axis of the system in the quadrupole 
approximation, 

c^{dh+) = Qj.^ - Qyy (23) 
c\dh^) = 2Q.,y. (24) 

In Fig. we show the corresponding gravity wave luminosity of the system, given by 

c^w = ^(gg'gg'). (25) 

We see that, as the inner NS cores merge, the gravity wave luminosity peaks for both runs, with the characteristic 
frequency of the waves increasing like (twice) the rotation frequency of the system. This frequency increase is more 
rapid in the PN case, since the inspiral is faster. 

After t ~ 30, the evolution of the N binary is rather straightforward. A triaxial object is formed at the center of the 
system, with spiral outflows emanating from the outer parts of each star. The spiral arms remain coherent for several 
windings before slowly dissipating, and finally leaving a low-density halo of material in the region r/R ~ 2 — 15. 
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During this time, the central triaxial object acts as the predominant source for the gravity waves as it spins down, 
leading to a characteristic damped oscillatory signature, at a luminosity approximately 1/30 that of the peak. The 
rise in central density from the initial value at i = to the final value at t = 80 is consistent with what is expected 
from the mass-radius relation for a Newtonian polytrope with F = 3. 

This simple picture, which is familiar from many previous Newtonian simulations, is seen to break down when IPN 
effects are taken into account. As is clear from Fig. || = 12), just prior to the final coalescence, IPN effects cause 
the long axis of each star to rotate forward relative to the binary axis, so that the inner part of each star leads the 
center of mass in the orbital rotation. This dynamical tidal lag is expected from the rapidly changing tidal forces 
during the final inspiral phase (LRS5). It is not to be confused with the tidal lag produced by viscous dissipation 
in nonsynchronized binaries (see, e.g., Ref. fs^). The dynamical tidal lag angle can be estimated analytically for a 
Newtonian binary whose orbit decays slowly by gravitational wave emission. Using Eq. (9.21) of LRS5, we estimate 
a lag angle at — 0.01 for 1/c^ = 0.16 and r ~ 2R. This is in agreement with the very small lag angle observed in our 
N run (barely visible at t — 12 in Fig. |l|). In contrast, from our FN run, we find at — 0.14, indicating that the more 
rapid inspiral can dramatically increase this effect. 

As the FN merger proceeds, material from the leading edge of each star wraps around the other, so that the cores 
simply slide past each other instead of striking more nearly head-on as in the N case. As this happens around i ~ 25, 
the maximum density drops slightly, and the gravity wave luminosity rises again, reaching a second peak at t ~ 35, 
with a maximum luminosity L2 = 0.65 Li compared to the first peak of luminosity Li. A cursory examination of Fig. |^ 
reveals a highly asymmetric, triaxial configuration near this time. The subsequent oscillations of the two cores in their 
sliding motion against each other damp out rather quickly, and the central object becomes more nearly axisymmetric 
while the maximum density rises again. A third peak of maximum luminosity L3 = 0.15 Li is clearly visible near 
f ~ 51, as is another very slight drop in the central density at that time, and a fourth, much smaller luminosity peak 
occurs at < ~ 72. 

To better understand this oscillation of the merger, and the corresponding modulation of the gravitational radiation 
waveforms, we show in Fig. ^ja comparison between the gravity wave luminosity and the ratio of the principal moments 
of inertia of the central object in the FN run. As can be seen clearly, the two quantities are strongly correlated. If 
we ignore the details of the internal motion of the fluid, it may be tempting to model the late-time behavior of the 
remnant in terms of a simple quadrupole (Z = 2 f-mode) oscillation of a rapidly and uniformly rotating single star. 
Adopting an average value for the angular velocity of the central object, 0.^ = 0.4, and using Eq. (3.30) of LRS5 for 
the frequency of the quadrupole oscillation of a compressible Maclaurin spheroid, we obtain a frequency a — 0.38, 
which gives us a modulation period Tmod = 16.6, very close to what we observe in Figs. 6 and 7. 

The occurrence of a second peak in the gravity wave luminosity can also be seen in the FN calculations presented 
in Ref. for polytropes with F — 2.6, but the second peak appears considerably less pronounced for F = 2.6 
than for F = 3. This may simply result from the higher central concentration of objects with lower values of F, 
which decreases the emission of gravitational radiation for a quadrupole deformation of given amplitude. Grid-based 
Newtonian calculations by Ruffert et al. |3l] for nonsynchronized binaries with a different EOS also show a second 
peak in the gravity- wave luminosity. For Newtonian systems with F ^ 2.2, the merger remnant evolves quickly to 
axisymmetry and the emission of gravitational radiation stops abruptly after the first peak (cf. RSI and RS2). 



C. The Final Merger Product 

In Fig. we show density contours of the central merger remnant in both the equatorial and vertical planes. For 
the N run, the remnant is shown a,t t — 80, which is at the end of the calculation. For the FN run, we show the 
remnant at t = 55, which corresponds to the third gravity wave luminosity peak, and at t = 80, the end of the 
simulation and close to a gravity wave luminosity minimum. Axes for the contour plots are aligned with the principal 
axes of the remnant. A summary of values for the principal axes and moments of inertia for the three configurations 
is presented in Table |. 

We see that the final remnant in the FN case is larger and more centrally condensed than in the N case, with a 
higher degree of fiattening in the vertical direction. This is in part because in the FN case less mass and angular 
momentum is extracted from the central region and deposited in the halo. Figures ^ and ^ show the evolution of 
the angular momentum of the various components in both runs. In the N case, most of the angular momentum lost 
by the remnant has gone into the halo. In the FN case, about equal amounts of angular momentum are lost to the 
halo and to the gravity waves. 

Nevertheless, the axis ratio a2/ai in the equatorial plane is approximately the same for the N run at t = 80 and 
for the FN run at t = 55 and at t = 80, indicating a reasonably constant shape for the outermost region. Further 
comparison between the N and FN remnants, however, shows that their interior structures are remarkably different. 
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In the PN remnant, the isodensity surfaces do not maintain a consistent orientation or shape as we move from the 
center to the equator of the remnant, indicating that the structure of the remnant is much more complex than that 
of a self-similar ellipsoid. Gravity-wave luminosity peaks are seen to occur when the inner and outer contours are 
aligned, leading to a larger net quadrupole moment (this is nearly the case at t = 55 in Fig. H). Minima occur when 
the orientations lie at right angles, as can be seen near t = 80 for the PN run in Fig. |^. 

In Fig. |l2|, we show the radial mass and rotational velocity profiles of the remnant. Horizontal cuts through the 
matter indicate that the rotation is cylindrical, with rotational velocity a function only of the distance from the 
rotation axis, independent of height relative to the equatorial plane (the same type of rotation profile has been 
obtained from strictly Newtonian calculations; see RSI). Neither case gives a rigid rotation law. The angular velocity 
of the N remnant shows a slight increase for r > 1.1 i?, whereas the PN run shows a decreasing angular velocity at the 
same point. Thus, both exhibit signs of differential rotation, but in opposite directions. We find that the centrifugal 
acceleration and gravitational acceleration become equal at the outer edge of the remnant for both cases, at r ~ 1.6 R 
and r ~ 1.85 R for the N and PN runs, respectively. This is in good agreement with the morphology of the remnants 
seen in Fig. ||, where a noticeable cusp-like deformation is visible in the outermost density contours near the equator 
in the vertical plane. We conclude that in both runs, the final remnant is maximally and differentially rotating. 

The rest mass of the N remnant at t = 80 is Mr = 1.73 M, while that of the PN remnant Mr = 1.90 M. The 
remaining mass, 0.27 M for the N run and 0.10 Af for the PN run, has been shed during the coalescence, forming 
the spiral arms seen in the middle panels of Figs. |^ and ^j. These spiral arms later merge to form a halo of matter 
around the central remnant. With a crude linear extrapolation from a halo mass of M^ — 0.27 M for the N run, with 
l/^iPAT = 0, and Mh = 0.10 M for the PN, run with l/c\pj^ = 0.05, we might expect that, for physically reasonable 
NS with 1/c^ ~ 0.15 — 0.20, the vast majority of the mass will remain in the central remnant. However, this result 
may be crucially dependent on our choice of initial spins and the EOS, and it is limited by the restrictions we have 
placed on the magnitude of the IPN corrections. It should also be noted that fully GR calculations of the coalescence 
of NS with a r = 2 EOS suggest that significant mass loss occurs even for extremely compact NS pO| . 



D. The Final Fate of the Remnant 



By their very nature our calculations cannot address directly the question of whether the NS merger remnant will 
collapse to form a BH. Indeed the parameters of our PN run were chosen so that all IPN quantities remain small 
throughout the evolution, which, for F ^ 4/3, guarantees stability. This can be verified directly by checking, for 
example, that the mass distribution in the final merger remnant remains everywhere well outside the corresponding 
horizon radius (see Fig. |l^). However, given some of the general properties of the merger remnant as determined by 
our calculations, we can ask whether an object with similar properties, but with a more realistic EOS and higher 
compactness, would still remain stable to collapse in full GR. For the coalescence of two 1.4 NS with realistic stiff 
EOS, it is by no means certain that the core of the final merged configuration will collapse on a dynamical timescale 



to form a BH (see Refs. |21 5q] for recent discussions 



The final fate of a NS binary merger in full GR depends not only on the NS EOS and compactness, but also on the 
rotational state of the merger remnant. It has been suggested, for example, that the Kerr parameter = Jr/Mg^ 
of the remnant may exceed unity for extremely stiff EOS |Q. This does not appear to be the case, at least for our 
choice of EOS. In Fig. we show the evolution of the Kerr parameter throughout the entire coalescence, including 
only particles for which the rest-mass density satisfies > 0.005. This cut includes essentially all matter in the initial 
stages, and effectively cuts out particles in the spiral outfiow once the coalescence begins, as well as those remaining 
in the outer halo at the end. We see that is very near unity just prior to the final merger, but, in contrast to 
what has been assumed in some previous studies |58), it decreases significantly during the final coalescence. The 
decrease occurs mainly during periods of maximum gravity-wave luminosity, as angular momentum is radiated away, 
and during the mass-shedding phase after t ~ 20, since angular momentum is transferred from the core to the outside 
spiral outfiow. By the end of the PN run, has decreased to ~ 0.7, well below unity, and certainly not large enough 
to prevent collapse. The final value of the Kerr parameter for the PN run, = 0.70, is considerably greater than 
that of the N run, = 0.47. The difference is attributable to the greater mass ejected in the N run, which carries off 
a significant fraction of the angular momentum of the system (see Figs. |l^ and pi] ). 

Quite apart from considerations of the Kerr parameter, the rapidly rotating core may be dynamically stable. Indeed, 
most stiff NS EOS (including the recent "AU" and "UU" EOS of Ref. ^) allow stable, maximally rotating NS with 
baryonic masses exceeding 3Mq i.e., well above the mass of the final merger core (which is 1.9 M ~ 2.85 M0 
for M = 1.5 Mq in our PN calculation; see Fig. |l2|). Differential rotation (not taken into account in the calculations 
of Ref. 1^) can further increase this maximum stable mass very significantly (see (5^). For slowly rotating stars, 
the same EOS give maximum stable baryonic masses in the range 2.5 — 3 Mq, implying that the core would probably 
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(but not certainly) collapse to a BH in the absence of rotational support. 

If the final merger remnant is being stabilized against collapse by rotation, one must then consider ways in which it 
may subsequently loose angular momentum. Further reduction of the angular momentum of the core by gravitational 
radiation or dynamical instabilities cannot occur, since, at the end of the dynamical coalescence, the core is, by 
definition, dynamically stable and nearly axisymmetric (i.e., no longer radiating gravity waves; see Figs. 6 and 7). 
The development of a secular bar-mode instability (a quadrupole mode growing unstably on the viscous dissipation 
timescale; see LRSl and LRS4) has been discussed as a way of reducing the angular momentum of a rapidly rotating 
compact object |6l|| . However, this cannot occur either for a binary merger remnant because, if the remnant were 
rotating fast enough to be secularly unstable, it would still be triaxial (Recall, for example, that the point of bifurcation 
of the classical Maclaurin spheroid sequence into the Jacobi ellipsoid sequence coincides with the onset of secular 
instability for Maclaurin spheroids; see, e.g., |Q and LRSl). Note that other processes, such as electromagnetic 
radiation or neutrino emission, which may also lead to angular momentum losses, take place on timescales much longer 
than the dynamical timescale (see, e.g., Ref. |6^, where it is shown that neutrino emission is probably negligible). 
These processes are therefore decoupled from the hydrodynamics of the coalescence. Unfortunately their study is 
plagued by many fundamental uncertainties in the microphysics. 

IV. SUMMARY AND DIRECTIONS FOR FUTURE WORK 

Using a Lagrangian, SPH-based adaptation of the BDS PN formalism for hydrodynamics, we have calculated the 
merger of a coalescing NS binary including IPN and gravitational radiation reaction effects. We have also developed 
a method for computing accurate, quasi-equilibrium initial data for coalescing binaries in PN gravity, improving upon 
previous calculations that used nonequilibrium initial conditions containing unperturbed, spherical stars. 

We have confirmed that PN corrections to gravity cause the binary inspiral to become dynamical at larger binary 
separation compared to what is predicted in the Newtonian limit. In calculations using Newtonian gravity, but 
including the effects of the gravitational radiation reaction, we have found that the inspiral rate just prior to merging 
agrees well with the predictions of semi-analytic models using compressible ellipsoids as trial functions in an energy 
variational method (LRS). 

Using a hybrid formalism where radiation reaction is treated realistically but IPN effects are reduced in amplitude 
so as to remain numerically tractable, we have compared the hydrodynamic coalescence of binary NS systems in 
Newtonian and PN gravity. We find that IPN effects lead to important qualitative differences in the hydrodynamic 
behavior and in the gravitational radiation waveforms and luminosities. In Newtonian gravity, the merger of two 
equal-mass, F = 3 polytropic NS produces a single peak in the gravity-wave luminosity, followed by an exponentially 
decaying signal. In PN gravity, we see a strong quadrupole oscillation of the remnant immediately after coalescence, 
which leads to several additional peaks in the gravity-wave luminosity. In both Newtonian and PN gravity, the final 
merger remnant is found to be maximally rotating and nearly axisymmetric. Even for realistic NS EOS and in full 
GR, this configuration is expected to be stable against gravitational collapse to a BH on a dynamical timescale. The 
amount of mass ejected into an outer halo by the rotational instability developing during the final merger decreases 
substantially when IPN effects are included, and we suggest that, for realistic NS models, essentially no mass might 
be ejected, so that the total baryonic mass of the system remains entirely in the central remnant (though this result 
is hardly a certainty). 

Our study is naturally beginning with PN calculations for equal-mass, corotating binaries with a simple polytropic 
EOS. This allows us to compare our results directly to previous Newtonian calculations performed with the same 
set of assumptions (RS, [ p6p7| , |3^ ). The dependence of our results on the NS EOS will be studied in future papers 
by varying the adiabatic exponent F (in the range F ~ 2 — 4 apphcable to NS; see, e.g., LRS3) and by performing 
additional runs with more realistic tabulated NS EOS. In particular, we will consider the Lattimer-Swesty EOS [ |63[ . 
This EOS includes high-temperature effects (which can be significant in the outermost, low-density regions of some 
NS mergers) and has also been employed in several previous Newtonian studies [p3[|, to which we want to compare 
our results. Even with the lowest available value of the nuclear compressibility {K = 180 MeV), the Lattimer-Swesty 
EOS is relatively stiff (effective F ~ 2.5 for a 1.4 Mq NS). The latest microscopic NS EOS, constrained by nucleon 
scattering data and the binding of light nuclei, and incorporating three-body forces, are even stiffer (effective F > 3; 
see, e.g., Ref. ||6^ for a summary, and Ref. for the latest version). We will use several of these recent EOS, in 
tabulated form, to perform additional, more realistic calculations. More schematic EOS based on exotic states of 
matter, such as pion condensates or strange quark matter, can be much softer (T 2 and maximum stable masses not 
much above 1.4 Af©). We will not consider such soft EOS in our calculations, since they render the PN approximation 
invalid. Note that several observations in progress may have already ruled them out (e.g., from the large measured 
mass of the NS in Vela X-1; 
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Using our PN SPH code we will also study the dependence of the hydrodynamics and gravitational wave emission 
on the binary mass ratio q. Neutron star masses derived from observations of binary radio pulsars are all consistent 
with a remarkably narrow underlying Gaussian mass distribution with Mg = 1.35±O.O4M0 |66j. The largest observed 
departure from g = 1 in any known binary pulsar with likely NS companion is currently q — 1.386/1.442 — 0.96 for the 
Hulse- Taylor pulsar PSR B1913+16 ||6^. Although the equal-mass case is clearly important, one should not conclude 
from these observations that it is unnecessary to consider coalescing NS binaries with unequal-mass components. 
Indeed, it cannot be excluded that other binary NS systems (that may not be observable as binary pulsars) could 
contain stars with significantly different masses. Moreover, Newtonian calculations of binary NS coalescence have 
shown that even very small departures from q = 1 can drastically affect the hydrodynamic evolution (RS2, [p9[). 

In future papers we will also perform PN SPH calculations for binary NS systems that are initially nonsynchronized . 
This is likely to be the case for real systems, since the tidal synchronization time in close NS binaries is probably 
always longer than the orbital decay time [^^. The methods of LRS can be used to construct approximate, quasi- 
equilibrium initial conditions for nonsynchronized coalescing binaries. For binaries that are far from synchronized, the 
final coalescence involves some new, complex hydrodynamic processes, and significant differences in the gravitational 
wave emission compared to the synchronized case, with an additional dependence of the gravitational radiation 



waveforms on the stellar spins P9 69 1. Moreover, the final fate of the merger may also be very different for initially 



nonsynchronized binaries, since the merger remnant may no longer be maximally rotating 



very gi 
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APPENDIX A: THE BDS 1 + 2.5 PN FORMALISM 

In the original Eulerian, PN formalism of Blanchet, Damour, and Schafer jiij, the key variables appearing in the 
hydrodynamic evolution equations are PN variants of the standard Newtonian quantities. Specifically, the coordinate 
rest-mass density and momentum per unit rest-mass Wi are given in terms of the proper rest-mass density p and 
the 4-velocity by 

r* = ^u°p (Al) 

(^1 + CUi, (A2) 

where h is the specific enthalpy of the fiuid. Assuming hereafter a polytropic equation of state, i.e. one for which the 
pressure is given by 

P*{r*) = krl, (A3) 

it is found that the specific enthalpy is given by 

h^k-^r^^-' = -^P^. (A4) 
I — 1 1 - 1 

It should be noted that is not the Newtonian pressure, but rather a IPN variant of it. 

The BDS formalism requires the solution of nine Poisson equations, one for the Newtonian gravitational potential 
?7,, seven for IPN corrections, and a final one to handle the gravitational radiation reaction. 

The equation for the gravitational potential is 

V^U, = -Ann. (A5) 

Note that with this sign convention, the gravitational potential is a positive quantity. The IPN correction potentials 
are given by 
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V^t/j = -Awr^w, (A6) 

V^Cj = -Airx'dsir^w,) (A7) 

= -^nrj. (A8) 

Note that the summation m Eq. ( |A6| ) runs over i ~ x,y,z, thus U2 7^ C/y in Eq. (|A8|). Using these, we define the 
quantity 

A^iU^ + ^Q-^x'dsUs. (A9) 



It is important to note that the volume integral of the source term of Eq. (A7) vanishes, assuming that t he o rigin 
is at the center of mass and momentum of the system, and thus it contains no monopole term. In Eq. ( |A§| ), the 
quantity 6 in the source term is one of three quantities which are assumed to be of order O(^). They are, assuming 



the equation of state Eq. (A3), and with w'^ — S^^WiW 



a = 2U,- T{-w^ + 3C/*) (AlO) 

The third derivative of the symmetric, trace-free (STF) quadrupole tensor, Q\^j is calculated from 

Pij ^2 J d^xr43w,djU^ - 2iy,^^ + x'w,d,jU^ - x'd^jUs] (A13) 

Q1? = lP^J + \Pj^ - \^^JPss. (AM) 

and is used in the source term for the radiation reaction potential C/5, of order O(^). This is calculated from the 
final Poisson equation, 

Uz = \G{R^Q^^h'djr,) (A15) 



V^R^-AirQf^x'djr^. (A16) 

Since we are dealing with the trace-free quadrupole tensor, it is easy to show that the volume integral of the source 
term of Eq. ( A16) also vanishes, for any mass distribution. 
Forces are defined by 

repress /, , '-'^ \ diP^ ^ P* ci / A 1 ■7\ 

Fl''^ = (^1 + d,u, + la,{/2 - ^wsd.As (A18) 

= \d,U5. (A19) 
Finally, the evolution system, in Eulerian form, is given by 



dtr, = d^ir^v') (A20) 
dtw^ = -v'dsw, + i^r'" + Fl'''' + J^r"', (A21) 



where the particle velocities are related to the specific coordinate momentum Wi by 



l-3H» + 72^^+5^^«C'- (A22) 



The quantities v and w will be referred to simply as the velocity and momentum vectors, respectively (see ||3T|). 
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In the SPH method, the evolution equations must be expressed in a Lagrangian form, given simply by 



X 

Wi 



V 

ppress 



F, 



IPN 



(A23) 
(A24) 



In BDS, there also appear evolution equations for the entropy and the pressure. The former, converted into 
Lagrangian form, states that entropy is a conserved quantity, and is handled in our code by the choice of AV, as in 
Eq. (^. The latter is not necessary here since we calculate the pressure directly from the density at each time step. 

Since the parameters a and /3, defined b y Eqs. (A10,All) become rather large for NS with ~ 0.05, we make 
some small adjustments to Eqs. ( A17 , A2^ ). We note that for a adiabatic exponent F > |, a is every wher e negative. 
To ensure that the pressure force always acts in the proper direction, we make a substitution in Eq. (A17), 



1 



- 1- 



"1 dip^ 



(A25) 



This new form is entirely equivalent to the one it replaces to IPN order, 
make the following substitution in Eq. ( A22 ), 



Wi 



Wi 



Similarly, /? is everywhere positive, so we 



(A26) 



APPENDIX B: RELAXATION METHODS 



1. PN Case 



Constructing hydrostatic equilibrium initial conditions in PN gravity is a much more difficult problem than in 
Newtonian gravity, primarily because of the complex relationship between the particle velocity and momentum. We 
get around this problem by implementing a multistage approach to the construction of relaxed configurations. 

First, we construct a series of hydrostatic equilibrium models for single F = 3 polytropes with increasing values of 
to gauge the effects of the PN corrections on the structure of the stars. Specifically, we construct relaxed models 
with compactness parameters of = 0.01 to 0.07, in steps of 0.01. 

In the relaxation procedure, spurious velocities arising from configurations adjusting toward equilibrium are ignored 
as sources for the force equations. Thus particles move during the relaxation, but the force exerted on each particle 
is th at of a static mass configuration. We thus solve all Poisson equations assuming w ^ 0, which eliminates 
Eqs. (A6,A7 ). In ad dition, the velocity terms in the definition of the IPN quantities a, (3, and d are removed from 
Eqs. (A1C,A11 A12). This greatly simplifies the equations giving us the set 



5 

repress 
plPN 



a ={2- 3F)[/, 



F - 1 7% 

3F - 2 
F- 1 



^ p* a 

— aa 



7^2 



1 - ^ I 



Wi 



repress ^ 



IPN 



Wi 
irela 



(Bl) 
(B2) 
(B3) 

(B4) 
(B5) 
(B6) 
(B7) 

(B8) 

(B9) 
(BIO) 
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where treiax is the relaxation time. 

To construct our first model, with \/c^ — 0.01 we start from a Newtonian F = 3 polytrope and let it relax to an 
equilibrium configuration. Then, using the maximum particle radius Rmax, we adjust the radial position, smoothing 
length, and specific entropy of all particles according to 

(Bll) 



^max 

hm - ^ (B12) 
-fir, 



34-3r 



km * ^mRmax ■ (-E^l^) 

Velocities are set to zero at the end of this rescaling. This new configuration is relaxed again, and the process 
is repeated until convergence is achieved. For the F = 3 models, we rescaled every t = 5.0, with a relaxation time 
treiax = 1-0, The final profile is used as the initial test configuration of the next model, which is then relaxed iteratively 
as described above. 

In addition to F = 3 models, we also computed a sequence of single FN polytropes with F — 5/3, and tested their 
stability. FN effects shouldmake the star unstable to gravitational collapse when 1/c^ ^ 0.141 for F = 5/3 [Q. We 
tested 1/c^ values in steps of 1/c^ = 0.02, until we reached 0.10, at which point we halved the step size until we 
reached 1/c^ = 0.13. To make the relaxation overdamped, we reduced the rescaling time to t = 2.0, with treiax ~ 1-0. 
It was found that 1/c^ = 0.13 is always unstable, collapsing inward uncontrollably, no matter how short the rescaling 
time. This agrees well with the theoretical prediction when we account for the magnitude of the IFN corrections we 



deal with, and the approximations made in the analytic treatment. In Fig. 14, we show the time evolution of the 
specific entropy k for both sequences, taken as a ratio with the Newtonian value of the specific entropy derived from 
the Lane-Emden equation. We see a gradual increase of k as the compactness is increased, in both sequences, until 
we get to 1/c^ = 0.12 for F — 5/3, for which k is 50% larger than the corresponding value for 1/c^ = 0.11. 

Parameters for the single star sequences are shown in Table ||. Radial profiles of the density, as well as all important 



IPN quantities are shown in Figs. 15 and |16|. We see in the F = 5/3 case that increasing the compactness increases the 
central concentration of the model, which can be seen in the factor of 2 increase in central density. For compactnesses 
near the stability limit, we see that a, f3, and S are all of order unity. A different behavior is seen in the F = 3 case, 
for which the internal structure of the star remains almost unchanged as the IPN order parameters get large. We 
see that a and f3 both get relatively large for more compact models, but 6 is rather small, since the potential and 
pressure terms cancel each other to some extent. 

A comparison of the mass profile for the F = 3 polytrope with — 0.05, the model used in the PN dynamical 
simulation, to a direct Runge-Kutta integration of the IPN structure equations in spherical symmetry is shown in 
Fig. ||. We see excellent agreement, except at the outer edge of the star, where surface effects alter the SPH mass 
profile slightly. This results from a layer of particles developing at the surface of the stars, with a slight density 
decrement immediately within, but involves only a very small fraction of the total mass of the system. Since our 
method restricts particle positions to r/i? < 1, we see that the density falls to zero slightly outside this point because 
of the finite size of the SPH smoothing kernel. 

Once these single star configurations were complete, the resulting stars were placed in duplicate in a binary config- 
uration, which was assumed to be in a state of synchronized rotation, i.e., the velocity of every SPH particle is given 
as a function of position by 

Wo = l1 X f. (B14) 

The main difficulty in relaxing PN configurations is in the interplay between v and w, which not only differ in 
magnitude but also in direction. Thus, one or the other can be relaxed in the corotating frame, but not both. Here 
V was assumed to be zero in the corotating frame for a relaxed configuration, satisfying the equation above. 

We created a method to calculate wo{vo), which is not invertible in closed form. As can be seen from Eq. ( [A2^ ), 



the relationship between particle velocity and momentum is a function of several potentials at the particle position, 



through the term containing Ai. Since Ai is itself a function of w (see Eq. A9), and vice versa, we need to solve 
consistently for both. It was found to be best to use an iterative procedure, which alternately solves for w and then 
uses these trial values in the source terms of the relevant Poisson equations. 

In the initial step, using known values of vq, we first approximate wq by the equations 

* = *(l + ^ + ^fl + M\ (B16) 
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The computed value of w enters into the source terms of both Ui and Ci (Eqs. A6 A7). Using these two potentials, 
we calculate Ai and P (Eqs. A£, All ), and recalculate a new approximation to wq, denoted Wnew, from the previous 
one, Wold, by an iterative method, using only ^ of the correction to avoid overshooting, thus 




It was found that, for the models we tested, about ten iterations would give c onve rgence to within 1 part in 10^ to the 
correct value of v when compared to the value of v{'Wnew) calculated by Eq. ( A22). For every timestep afterwards, we 
followed the same iteration procedure, and about six iterations were found to produce the same convergence to the 
proper values. 

Once convergence to an acceptable solution was found, forces were calculated, and ii was estimated by finite 
differencing, 

v{w{t + dt))-v{w{t)) 

Vforce = • l^J^oj 

We relax the binary models at fixed center-of-mass separation r, in the corotating frame, adjusting such that the 
inward force of gravity is balanced exactly by the centrifugal force. At every time step, we calculate 



n = d ^ " , (B19) 

where Fin refers to the net inward force on each component of the binary. Particle velocities are advanced according 
to 

V = Vforce- +^^r. (B20) 

''relax 

After every time step, the two stars were adjusted slightly to maintain a center of mass separation at the desired 
value. 



2. Newtonian Case 



In the regime where the dynamical timescale of the neutron stars is much smaller than the characteristic timescale 
for gravitational radiation, we expect the stars to evolve through a series of quasi-equilibrium configurations. If 
synchronized rotation is assumed, these equilibrium configurations can be constructed by adding a centrifugal force 
and drag term to the acceleration equation, giving us 



where the centrifugal potential is given by 



^;''^"^™-v,($ + $™0-T^, (B21) 

'-'relax 

'^rot^^n^x^ +y^). (B22) 

The relaxation timescale, treiax is set initially to 1.0, close to the value required for critical damping of oscillations 
(RSI). For the purposes of relaxation, AV and the radiation back-reaction, which are both time-asymmetric, are 
ignored. In addition, during the relaxation, we ignore the distinction between velocity and momentum vectors in 



Eq. ([Tsl), taking v = w. The rate of rotation is calculated as in the PN case by Eq. (B19). Once the binary has 
relaxed to a suitable initial configuration, it is set in motion, and we commence the dynamical run. Initial velocities 
are given by 

Wx — ^^y, Wy — (B23) 



and V is calculated from w by Eq. (A22). In the point mass limit, this would reduce to Eq. (35) of Ruffcrt, Janka, 
and Schafer who use 

16 M3 , ^ 

vr^--— B24 
5 r'' 

as their initial condition. 
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FIG. 1. Evolution of the system in the N run. Projections of a random subset of 20% of all SPH particles onto the orbital 
(x-y) plane are shown at various times. The orbital motion is counter-clockwise. Units are such that G = M = R = 1, where 

M and R are the mass and radius of a single, spherical NS. Note that the development of a mass-shedding instability after 
t ~ 25, and the rapid contraction of the remnant toward an axisymmetric state at late times. 
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FIG. 3. Evolution of the binary center of mass separation during the inspiral phase for the two calculations. The solid line 
is for the PN run, the dashed line for the N run. The horizontal line represents the dynamical stability limit for a Newtonian, 
equilibrium binary, at r ~ 2.7 R. It appears as a break in the inspiral rate of the Newtonian binary, whereas the PN binary 
inspiral becomes dynamical at significantly greater separation. 
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FIG. 4. Same as Fig. but focusing on the early inspiral of the Newtonian binary. The sohd hne is the result from the SPH 
calculation (N run). The dashed line shows the point-mass approximation, the dash-dotted and dotted lines the approximations 
for two spheres and two ellipsoids, respectively. See text for details. The point mass approximation clearly fails when tidal 
interactions become significant, but note the excellent agreement with semi-analytic results for extended stars before the ISCO 
is encountered. 
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FIG. 5. Evolution of the maximum density in the two coalescence calculations. The upper curve is for the PN run, the lower 
curve for the N run. The sharp decline in density at t ~ 15 occurs as the two NS are tidally disrupted, followed by a larger 
increase as they coalesce. 
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FIG. 6. Gravity wave signatures for the two coalescence runs. The waveforms are calculated for an observer at a distance 
d along the rotation axis. The solid line shows the /i+ polarization, the dashed line the hx polarization (see Eqs. p3|j24| ). At 
late times in the N run the waveforms show a simple, exponentially damped oscillation, whereas in the PN run an additional 
large-amplitude modulation is apparent. 
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FIG. 7. Gravity wave luminosity for the two coalescence runs (see Eq. |2^). The solid line is for the PN run, the dashed line 
for the N run. The peak luminosity in the PN run is smaller than that of the N run, but secondary peaks occur at t ~ 35, 50, 
and 70. 
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FIG. 8. The ratio of the principal moments of inertia in the equatorial plane for the PN merger remnant, compared to the 
gravity wave luminosity at late times. The times of maximum elongation correspond to maxima in the gravity wave luminosity, 
and to decreases in the maximum density in Fig. H at t ~ 35 and t ~ 50 (and less clearly at t ~ 70). 
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FIG. 9. Density contours of the merger remnants. The top frames show the PN remnant at t = 55, the middle ones show the 
same remnant at t = 80, and the lower ones show the N remnant at t = 80. The left frames show a cut through the equatorial 
plane, the right frames through the vertical plane (containing the rotation axis). Contours arc logarithmic, ten per decade, 
starting from the maximum density of (r»)max = 0.567 for the PN run at t = 55, {rt)max = 0.608 for the PN run at t = 80, 
and {rt.)max = 0.518 for the N run at t = 80. The axes have been rotated to fall along the principal axes of the remnant. Note 
the cusp-like shape of the contours near the equator in the vertical plane, indicating maximal rotation. 
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FIG. 10. Evolution of the angular momentum in various components in the N run. Here Jtot is the total angular momentum 
in the system, Jr is for the inner remnant (defined by the condition r* > 0.005, which includes the entire binary initially, but 
only the inner remnant at later times), and Jh is for the outer halo (so that Jtot = Jr + Jh)- The dotted line shows the initial 
angular momentum of the system. 
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FIG. 11. Evolution of the angular momentum in various components in the PN run. Conventions are as in Fig. |l^. 
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FIG. 12. Enclosed rest mciss and radially averaged rotational velocity profiles of the final merger remnant at t=80 for the 
two runs. Here, Vcyi is the distance from the rotation axis, while r is the radius from center. In both plots, the solid line is for 
the PN run, the dotted line for the N run. The dashed line shows the radius for a Kerr black hole with a = 0.7 (the value we 
find for the PN run at t=80). For the rotational profile, we show only the data for —0.1 < a < 0.1, all other horizontal cuts 
yielding similar profiles extending to smaller radii. 
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FIG. 13. The evolution ol the Kerr parameter = cJr/Mg^, for the inner remnant in the PN run (sohd hue) and in the N 
run (dashed hue). At no time do we have > 1- The inner remnant (or core) is defined by the same density cut as in Fig. |l^. 
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FIG. 14. Results of SPH relaxation calculations for single stars. The ratio of the PN specific entropy fc to the Newtonian 
value fejv, is shown for both P = 5/3 and P = 3, computed for sequences of increasing compactness The dotted lines 

give the final value for each case, which was used as the initial value for the next relaxation. For P = 5/3, we see that for 
Ijc? > 0.12, fe/fcjv increases without bounds, indicating instability. For P = 3 and 1/c^ > 0.07, the IPN approximation breaks 
down. 
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TABLE I. Properties of the Merger Remnants. Units are such that G = M = R = 1, where M and R are the mass and 
radius of a single, spherical NS. Here, Mr is the rest mass of the remnant, Mgr is its gravitational mass, Jr is its total angular 
momentum, flc and fleq arc the angular rotation velocities at the center and at the equator, and the Oi's and li's are the 
principle axes and moments of inertia. 
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TABLE II. Parameters for Single Star Models. For each model, we list the compactness parameter 1/c^, the ratio of the PN 
specific entropy k to the Newtonian value fcjv, and the central values of density r*, and the dimensionless ratios P/r*(? and 



lie- 


k/kN 




{P/r.c% 


{U./c% 


r = 5/3 


0.02 


1.177 


1.201 


0.0111 


0.0438 


0.04 


1.281 


1.292 


0.0257 


0.0893 


0.06 


1.421 


1.452 


0.0464 


0.1377 


0.08 


1.634 


1.708 


0.0792 


0.1902 


0.10 


1.879 


1.976 


0.1255 


0.2470 


0.11 


2.198 


2.336 


0.1806 


0.2820 


0.12 


3.400 


2.295 


0.3011 


0.2968 


r = 3 


0.01 


1.403 


0.3822 


0.0052 


0.0165 


0.02 


1.553 


0.3818 


0.0114 


0.0326 


0.03 


1.649 


0.3882 


0.0187 


0.0486 


0.04 


1.780 


0.3948 


0.0280 


0.0649 


0.05 


1.918 


0.4051 


0.0397 


0.0813 


0.06 


2.084 


0.4170 


0.0549 


0.0989 


0.07 


2.262 


0.4321 


0.0746 


0.1154 
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